package triangulation;

import java.awt.Point;
import java.nio.FloatBuffer;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.Random;

import org.lwjgl.BufferUtils;

import computation.Curvatures;
import computation.Geometry;
import computation.Smoothing;

import color.CLRConst;
import color.CLRConvert;

import delaunay_triangulation.Delaunay_Triangulation;
import delaunay_triangulation.Point_dt;
import delaunay_triangulation.Triangle_dt;

import primatives.PointND;
import primatives.Vector;
import quadtree.QuadBoundingBox;
import quadtree.QuadTree;
import ui.UIConst;
import util.GeoEdge;
import util.GeoTriangle;
import util.GeoVertex;
import util.Matrix_3X3;
import util.STLTriangle;

import static org.lwjgl.opengl.GL11.*;
import static org.lwjgl.opengl.GL15.*;

public class SimpleTriangulation implements Triangulation
{
	/* TODO:
	 * - experiment with display lists recompling and creating new VBOs for changing things like color and geometry.
	 * - look into ARBVertexBufferObject.glBindBufferARB() for modifying the buffers directly rather than regenerating each time
	 * - finish implementing Taubin smoothing
	 */

	private static double EPSILON = 0.000001; // extremely small non-zero value for functions that need a positive number
	
	private static final int TRIANGLE_SIZE = 3; // used to tell the VBO the number of vertices per face
	private static final int VERTEX_SIZE = 3; // used to tell the VBO the number of floats per vertex
	private static final int EDGE_SIZE = 2; // used to tell the VBO the number of vertices per edge
	
	protected QuadTree _tree;
	protected Delaunay_Triangulation _delaunay;
	public int CURRENT_SMOOTHNESS_LEVEL;
	public int MAX_SMOOTHNESS_LEVEL = 0;

	public DelaunayMapping mapping;
	protected int _numOfSelectedTriangles;
	
	protected boolean _lightsEnabled;
	protected boolean _wireframe;
	protected boolean _point;
	protected boolean _renderToRender;
	protected SELECTION_MODE _mode;
	
	protected double roughness_std;
	protected double roughness_avg;
	protected double roughness_max;
	
	// color range variables, ensure that these are always positive numbers
	private double[] color_scale_range;
	private double[] data_scale_range;
	private double data_offset; // offset value that shifts the actual data range into the color range
	public SCALING_MODE colorScaleMode;
	
	// window querying points, these points represent the actual app window
	public Point topleft;
	public Point bottomright;
	
	public boolean isDragging;
	public boolean allowDragging;
	public boolean updateComplete;
	public boolean updateGeometry; 
	public boolean updateColors;
	public boolean updateNormalLines;
	
	public boolean perVertexStrikeDip = true;
	public boolean perTriangleStrikeDip = false;
	public boolean showNormals = true; // flag for displaying the normals
	
	// TODO: do not use these for now, figure out how to properly do per vertex lighting and then continue
	public boolean usePerTriangleNormal = false;
	public boolean usePerVertexNormal = false;
	
	public boolean bStrikeVersion1 = false; // flag for switching between the strike and dip versions
	
	// strike and dip of the laser head, ignore these parameters if they are both negative
	public double dLaserStrike = -1.0;
	public double dLaserDip = -1.0;
	
	
	private float MESH_ZOOM;

	public enum SELECTION_MODE
	{
		NONE,
		STRIKE_ONLY,
		DIP_ONLY,
		STRIKE_DIP,
		STRIKE_DIP_SEARCH, // only this makes use of the triangle "selected" flag
		ROUGHNESS,
		PEAK_VALLEY_SADDLE,
		RANGE_QUERY,
		ROUGHNESS_MAP,
		CURVATURE_ANISOTROPIC,
		CURVATURE_ANISOTROPIC_AVG,
		CURVATURE_MEAN,
		CURVATURE_GAUSSIAN,
		CURVATURE_RING_MAPS,
		CYLINDER_CLUSTERING
	}
	
	public enum SCALING_MODE
	{
		NONE,
		LOG,
		SQRT
	}
	
	public SimpleTriangulation()  { init(); }
	
	public SimpleTriangulation(ArrayList<PointND> pointCloud) { init(); triangulate(pointCloud); }
	
	private void init()
	{
		mapping = null;
		
		// flags for dragging the mouse
		isDragging = false;
		allowDragging = true;
		
		updateComplete = false; // flag which triggers a complete regeneration of the geometry & painting
		updateGeometry = false; // flag which triggers a complete regeneration of the geometry 
		updateColors = false; // flag which triggers a repaint_triangulation
		updateNormalLines = false; // flag which triggers a regeneration of the normal lines
		
		// TODO: do not use these for now, figure out how to properly do per vertex lighting and then continue
		usePerTriangleNormal = false;
		usePerVertexNormal = false;

		_lightsEnabled = false;
		_wireframe = false;
		_point = false;
		_renderToRender = false; // flag to indicate if all our buffers+geometry is created and ready to be displayed
		_numOfSelectedTriangles = 0;
		_mode = SELECTION_MODE.NONE;
		colorScaleMode = SCALING_MODE.NONE;
		
		MESH_ZOOM = 1.0f; // a zooming variable does not actually modify the mesh
		CURRENT_SMOOTHNESS_LEVEL = 0; // controls the smoothness level of the mesh
		
		color_scale_range = new double[2];
		data_scale_range = new double[2];
		data_offset = 0.0;
	}

	// useful for communicating between triangulation and GeoCanvas what range the color scale is representing
	public void setColorScaleRange(double min, double max)
	{
		setDataRange(min, max);
		
		double offset = 0.0;
		if (min < 0)
		{
			offset = min*-1;
		}
		
		data_offset = offset;
		
		double clr_min = min + offset + EPSILON;
		double clr_max = max + offset + EPSILON;
		
		color_scale_range[0] = clr_min;
		color_scale_range[1] = clr_max;
	}
	public double[] getColorScaleRange()
	{
		double[] scaled_range = {color_scale_range[0], color_scale_range[1]};
		
		// color range does not account for log or sqrt scaling so scale them now
		if (isScaleMode_Log())
		{
			scaled_range[0] = Math.log(scaled_range[0]);
			scaled_range[1] = Math.log(scaled_range[1]);
		}
		else if (isScaleMode_Sqrt())
		{
			scaled_range[0] = Math.sqrt(scaled_range[0]);
			scaled_range[1] = Math.sqrt(scaled_range[1]);
		}
		
		return scaled_range;
	}
	
	public void setDataRange(double min, double max)
	{
		data_scale_range[0] = min;
		data_scale_range[1] = max;
	}
	public double[] getDataRange() { return data_scale_range; }
	
	public boolean showColorScale()
	{
		if (_mode == SELECTION_MODE.ROUGHNESS_MAP || _mode == SELECTION_MODE.CURVATURE_ANISOTROPIC || 
			_mode == SELECTION_MODE.CURVATURE_ANISOTROPIC_AVG || _mode == SELECTION_MODE.CURVATURE_MEAN || 
			_mode == SELECTION_MODE.CURVATURE_GAUSSIAN || _mode == SELECTION_MODE.CURVATURE_RING_MAPS) return true;
		boolean bScaleMode = _mode == SELECTION_MODE.STRIKE_ONLY || _mode == SELECTION_MODE.DIP_ONLY || _mode == SELECTION_MODE.STRIKE_DIP; 
		return bScaleMode || _numOfSelectedTriangles > 0;
	}
	
	// update the color scales according to the new mode
	public void resetColorScale_Curvatures()
	{
		// find the max and min
		double data_max = Double.MIN_VALUE;
		double data_min = Double.MAX_VALUE;
		for (GeoVertex g: mapping.vertices)
		{
			GeoVertex gv = mapping.getVertex(g.id, CURRENT_SMOOTHNESS_LEVEL);
			double value = 0.0;
			
			if (isMode_RoughnessMaps()) value = gv.curvatureRoughness;
			else if (isMode_CurvatureAnisotropic()) value = gv.curvatureAnisotropic;
			else if (isMode_CurvatureAnisotropicAvg()) value = gv.curvatureAnisotropicAveraged;
			else if (isMode_CurvatureMean()) value = gv.curvatureMean;
			else if (isMode_CurvatureGaussian()) value = gv.curvatureGaussian;
			else if (isMode_CurvatureRingMap()) value = gv.curvatureAvgRing;

			if (Double.isNaN(value)) value = 0.0; // TODO: fix all values that give a NaN!
			data_max = Math.max(value, data_max);
			data_min = Math.min(value, data_min);
		}
		setColorScaleRange(data_min, data_max);
	}
	
	private float _xRot = 0.0f;
	private float _yRot = 0.0f;
	private float _zRot = 0.0f;
	public void rotate(float x, float y, float z)
	{
		_xRot += x;
		_yRot += y;
		_zRot += z;
	}
	
	public void reset_rotations()
	{
		_xRot = 0.0f;
		_yRot = 0.0f;
		_zRot = 0.0f;
	}
	
	public double xBoundMin, xBoundMax, yBoundMin, yBoundMax, zBoundMin, zBoundMax;
	public double getBoundX_Legnth() { return xBoundMax - xBoundMin; }
	public double getBoundY_Legnth() { return yBoundMax - yBoundMin; }
	public double getBoundZ_Legnth() { return zBoundMax - zBoundMin; }

	public void triangulate(ArrayList<PointND> pointCloud)
	{
		// quad tree is really only needed if one requires selecting subsets of vertices for range querying
		if (UIConst.ENABLE_MOUSE_MESH_SELECTION)
		{
			System.out.print("Creating QuadTree...");
			_tree = new QuadTree(pointCloud);
			System.out.println("   done!");			
		}
		
		if (CURRENT_SMOOTHNESS_LEVEL > 0) smoothMesh_Reset();

		xBoundMin = yBoundMin = zBoundMin = Double.MAX_VALUE;
		xBoundMax = yBoundMax = zBoundMax = Double.MIN_VALUE;
		
		System.out.print("Triangulating...");
		Point_dt[] points = new Point_dt[pointCloud.size()];
		for (int i = 0; i < points.length; i++)
		{
			double x = pointCloud.get(i).get_coord(0);
			double y = pointCloud.get(i).get_coord(1);
			double z = pointCloud.get(i).get_coord(2);
			
			if (x > xBoundMax) xBoundMax = x;
			if (x < xBoundMin) xBoundMin = x;
			if (y > yBoundMax) yBoundMax = y;
			if (y < yBoundMin) yBoundMin = y;
			if (z > zBoundMax) zBoundMax = z;
			if (z < zBoundMin) zBoundMin = z;
			
			points[i] = new Point_dt(x, y, z);
		}

		_delaunay =  new Delaunay_Triangulation(points);
		System.out.println("   done!");
		
		System.out.println("Point Cloud Bounds: X(" + xBoundMin+","+ xBoundMax +")\tY("+yBoundMin+","+yBoundMax+")"+")\tZ("+zBoundMin+","+zBoundMax+")");

		_wireframe = false;
		_point = false;
		_renderToRender = false;
		
		_mode = SELECTION_MODE.NONE;
		_numOfSelectedTriangles = 0;
		
		createMapping();
	}
	
	public void createMapping()
	{
		System.out.print("Creating initial geometry mappings... ");
		if (mapping == null) mapping = new DelaunayMapping();
		mapping.clearMapping();
		mapping.createMapping(_delaunay);
		System.out.println("done!");
	}
	
	public Iterator<Triangle_dt> getTriangles() { return _delaunay.trianglesIterator(); }
	public Iterator<Point_dt> getVertices() { return _delaunay.verticesIterator(); }
	
	public int getNumberOfSelectedTriangles() { return _numOfSelectedTriangles; }
	public int getSize() { return _delaunay.size(); }
	public boolean isLoaded() { return _delaunay != null && _delaunay.size() > 0; }
	
	public void setMode(SELECTION_MODE m) { _mode = m; }
	public boolean isMode_StrikeOnly() { return _mode == SELECTION_MODE.STRIKE_ONLY; }
	public boolean isMode_DipOnly() { return _mode == SELECTION_MODE.DIP_ONLY; }
	public boolean isMode_StrikeDipOnly() { return _mode == SELECTION_MODE.STRIKE_DIP; }
	public boolean isMode_StrikeDipSearch() { return _mode == SELECTION_MODE.STRIKE_DIP_SEARCH; }
	
	public boolean isMode_Roughness() { return _mode == SELECTION_MODE.ROUGHNESS; }
	public boolean isMode_PeakValleySaddle() { return _mode == SELECTION_MODE.PEAK_VALLEY_SADDLE; }
	public boolean isMode_RangeQuery() { return _mode == SELECTION_MODE.RANGE_QUERY; }
	public boolean isMode_RoughnessMaps() { return _mode == SELECTION_MODE.ROUGHNESS_MAP; }
	public boolean isMode_CurvatureAnisotropic() { return _mode == SELECTION_MODE.CURVATURE_ANISOTROPIC; }
	public boolean isMode_CurvatureAnisotropicAvg() { return _mode == SELECTION_MODE.CURVATURE_ANISOTROPIC_AVG; }
	public boolean isMode_CurvatureMean() { return _mode == SELECTION_MODE.CURVATURE_MEAN; }
	public boolean isMode_CurvatureGaussian() { return _mode == SELECTION_MODE.CURVATURE_GAUSSIAN; }
	public boolean isMode_CurvatureRingMap() { return _mode == SELECTION_MODE.CURVATURE_RING_MAPS; }
	
	public boolean isScaleMode_None() { return colorScaleMode == SCALING_MODE.NONE; }
	public boolean isScaleMode_Log() { return colorScaleMode == SCALING_MODE.LOG; }
	public boolean isScaleMode_Sqrt() { return colorScaleMode == SCALING_MODE.SQRT; }
	
	public void setSmoothnessLevel(int level)
	{
		CURRENT_SMOOTHNESS_LEVEL = level;
		updateComplete = true;
	}
	public int getSmoothnessLevel() { return MAX_SMOOTHNESS_LEVEL; }
	public void setMeshScale(float s) { MESH_ZOOM = s; }
	public void zoom(float s) { MESH_ZOOM += s; }
	
	// lights control
	public boolean isLightsEnabled() { return _lightsEnabled; }
	public void enableLights() { _lightsEnabled = true; }
	public void disableLights() { _lightsEnabled = false; }

	// render types control
	public void setRenderFormat_Wireframe()
	{
		_wireframe = true;
		_point = false;
	}
	public void setRenderFormat_Points()
	{
		_wireframe = false;
		_point = true;
	}
	public void setRenderFormat_Polygon()
	{
		_wireframe = false;
		_point = false;
	}
	
	// computes normals per vertex, needs to called after createGeometry as we need the hash map of triangles first
	// averages the normals from the triangles which surround a given vertex
	public void computeVertexNormals()
	{
		System.out.print("Computing normals per vertex... ");
		for (GeoVertex gv: mapping.vertices)
		{
			float[] vertex_normal = new float[3];
			vertex_normal[0] = vertex_normal[1] = vertex_normal[2] = 0.0f;
			HashSet<Integer> triangkeKeys = gv.trianglekeys;
			for (Integer tKey: triangkeKeys)
			{
				GeoTriangle g = mapping.triangles[tKey];
				float[] n = g.getNormal();
				vertex_normal[0] += n[0];
				vertex_normal[1] += n[1];
				vertex_normal[2] += n[2];
			}
			
			vertex_normal[0] /= triangkeKeys.size();
			vertex_normal[1] /= triangkeKeys.size();
			vertex_normal[2] /= triangkeKeys.size();
			
			gv.setNormal(vertex_normal);
		}
		System.out.println(" done!");
	}

	private int vboHandle_Vertices;
	private int vboHandle_Normals;
	private int vboHandle_NormalLines;
	private int vboHandle_Colors;
	private FloatBuffer buf_vertices;
	private FloatBuffer buf_normals;
	private FloatBuffer buf_normal_lines; // a buffer of lines representing all the normals per triangle
	private FloatBuffer buf_colors;
	// called once when loading a completely new geometry, needed for VBO rendering
	// returns true if geometry has been created
	public boolean createGeometry()
	{
		if (this.isLoaded() && !_renderToRender && mapping.bCreated) // done once after the DT is completed
		{
			System.out.print("Creating initial OpenGL geometry ... ");
			final int numberTriangles = mapping.triangles.length;
			
			// create the buffers and maps
			buf_vertices = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 vertex per 1 triangle
			buf_normals = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 normal per 1 triangle
			buf_colors = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 color per 1 vertex
			buf_normal_lines = BufferUtils.createFloatBuffer(numberTriangles * VERTEX_SIZE * 2); // 3 floats per normal per triangle x2 as we want a line

			int triangleIndex = 0;
	    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
			while (iter.hasNext())
			{
				Triangle_dt dt = iter.next();
				Point_dt v0 = dt.p1();
				Point_dt v1 = dt.p2();
				Point_dt v2 = dt.p3();
				if (v0 != null && v1 != null && v2 != null)
				{
					// fill the vertex buffer
					float[] triangle = {(float)v0.x(), (float)v0.y(), (float)v0.z(),
								(float)v1.x(), (float)v1.y(), (float)v1.z(),
								(float)v2.x(), (float)v2.y(), (float)v2.z()};
					buf_vertices.put(triangle);
					
					// fill the normals buffer
					GeoTriangle gt = mapping.triangles[triangleIndex];
					GeoVertex gv0 = mapping.vertices[gt.vertexkeys[0]];
					GeoVertex gv1 = mapping.vertices[gt.vertexkeys[1]];
					GeoVertex gv2 = mapping.vertices[gt.vertexkeys[2]];
					
					buf_normals.put(Geometry.computeVertexNormal(gv0, mapping)); // per vertex normals yo
					buf_normals.put(Geometry.computeVertexNormal(gv1, mapping));
					buf_normals.put(Geometry.computeVertexNormal(gv2, mapping));
					
					float[] normal = Geometry.computeSurfaceNormal(gv0, gv1, gv2);
					gt.setNormal(normal);
					
					// compute the normal lines
					float[] triangleCentre = new float[3];
					triangleCentre[0] = (float)((v0.x() + v1.x() + v2.x())/3);
					triangleCentre[1] = (float)((v0.y() + v1.y() + v2.y())/3);
					triangleCentre[2] = (float)((v0.z() + v1.z() + v2.z())/3);
					buf_normal_lines.put(triangleCentre);
					float[] normal_line = Geometry.vectorPlus(triangleCentre, Geometry.vectorNormalize(normal));
					buf_normal_lines.put(normal_line);
					
//					double[] strike_n_dip;
//					if (bStrikeVersion1) strike_n_dip = computeStrikeNDipAngle(normal);
//					else strike_n_dip = computeStrikeNDipAngle2(normal);
					
					double[] strike_n_dip = Geometry.computeStrikeNDipAngle2(normal, dLaserStrike, dLaserDip);
					
					if (!Double.isNaN(strike_n_dip[1]) && !Double.isNaN(strike_n_dip[0]))
					{
						gt.StrikeDeg = strike_n_dip[0];
						gt.DipDeg = strike_n_dip[1];
						
//						System.out.println("["+normal[0] + "," + normal[1]+","+normal[2]+"]\t("+gt.StrikeDeg + "," + gt.DipDeg+")");
					}
					
					// set the intial color to be green
					buf_colors.put(CLRConst.GREEN);
					buf_colors.put(CLRConst.GREEN);
					buf_colors.put(CLRConst.GREEN);
					
					triangleIndex++;
				}
			}
			
			buf_vertices.flip();
			buf_normals.flip();
			buf_colors.flip();
			buf_normal_lines.flip();

			vboHandle_Vertices = glGenBuffers();
			vboHandle_Normals = glGenBuffers();
			vboHandle_Colors = glGenBuffers();
			vboHandle_NormalLines = glGenBuffers(); // normal lines
			
	        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Vertices);
	        glBufferData(GL_ARRAY_BUFFER, buf_vertices, GL_STATIC_DRAW);
	        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Normals);
	        glBufferData(GL_ARRAY_BUFFER, buf_normals, GL_STATIC_DRAW);
	        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Colors);
	        glBufferData(GL_ARRAY_BUFFER, buf_colors, GL_STATIC_DRAW);
	        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_NormalLines);
	        glBufferData(GL_ARRAY_BUFFER, buf_normal_lines, GL_STATIC_DRAW);
	        
	        glBindBuffer(GL_ARRAY_BUFFER, 0);
	        
	        
	        // compute any extra information here
	        computeRoughness(); // a measure of roughness via the dot product
	        computePeakValleySaddle(); // a measure of peak ratio
	        
			_renderToRender = true;
			System.out.println(" done!");
			
			System.out.println("Number of Vertices: " + _delaunay.size());
			System.out.println("Number of Triangles: " + _delaunay.trianglesSize());
			return true;
		}
		return false;
	}
	
	// completely regenerate everything
	public void regenerate_complete()
	{
		regenerate_geometry();
		regenerate_normal_lines();
		repaint_triangulation();
		System.out.println(); // extra line for good looks
		System.gc(); // clean up old memory
	}
	
	// regenerates all the geometry: vertices, faces and normals
	public void regenerate_geometry()
	{
		System.out.print("Regenerating geometry...");
		
		final int numberTriangles = mapping.triangles.length;
		_renderToRender = false;

		// ensure that we de-reference these buffers
		buf_vertices = null;
		buf_normals = null;
		
		// create the buffers and maps
		buf_vertices = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 vertex per 1 triangle
		buf_normals = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 normal per 1 triangle
		
		int triangleIndex = 0;
    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		while (iter.hasNext())
		{
			Triangle_dt dt = iter.next();
			Point_dt v0 = dt.p1();
			Point_dt v1 = dt.p2();
			Point_dt v2 = dt.p3();
			if (v0 != null && v1 != null && v2 != null)
			{
				// first get the geo triangle representation
				GeoTriangle gt = mapping.triangles[triangleIndex];

				// fill the vertex buffer with the correct vertices
				GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
				GeoVertex gv1 = mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
				GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);

				float[] triangle = {(float)gv0.x(), (float)gv0.y(), (float)gv0.z(),
						(float)gv1.x(), (float)gv1.y(), (float)gv1.z(),
						(float)gv2.x(), (float)gv2.y(), (float)gv2.z()};
				
				buf_vertices.put(triangle);
				
				// fill the normals buffer (we add the normal 3 times as OpenGL does lighting per vertex)
				buf_normals.put(Geometry.computeVertexNormal(gv0, mapping));
				buf_normals.put(Geometry.computeVertexNormal(gv1, mapping));
				buf_normals.put(Geometry.computeVertexNormal(gv2, mapping));
				

				
				float[] normal = computeSurfaceNormal(gv0, gv1, gv2); // find the normal = (x, y, z) = (East, North, Up)
				gt.setNormal(normal); // reset the normal so all other computations go off the new normals
//				buf_normals.put(normal);
//				buf_normals.put(normal);
//				buf_normals.put(normal);
				
				triangleIndex++;
			}
		}
		
		buf_vertices.flip();
		buf_normals.flip();

        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Vertices);
        glBufferData(GL_ARRAY_BUFFER, buf_vertices, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Normals);
        glBufferData(GL_ARRAY_BUFFER, buf_normals, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, 0);

		_renderToRender = true;
		System.out.println(" done!");
	}
	
	// regenerate the normal array
	public void regenerate_normals()
	{
		System.out.print("Regenerating normals...");
		
		final int numberTriangles = mapping.triangles.length;
		_renderToRender = false;
		
		// ensure that we de-reference these buffers
		buf_normals = null;
		
		// create the buffers and maps
		buf_normals = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 normal per 1 triangle
		
		int triangleIndex = 0;
    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		while (iter.hasNext())
		{
			Triangle_dt dt = iter.next();
			Point_dt v0 = dt.p1();
			Point_dt v1 = dt.p2();
			Point_dt v2 = dt.p3();
			if (v0 != null && v1 != null && v2 != null)
			{ 
				// first get all the key representations of the vertices
				GeoTriangle gt = mapping.triangles[triangleIndex];

				// fill the vertex buffer with the correct vertices
				GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
				GeoVertex gv1 = mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
				GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
				
				if (usePerTriangleNormal)
				{
					// fill the normals buffer (we add the normal 3 times as OpenGL does lighting per vertex)
					float[] normal = computeSurfaceNormal(gv0, gv1, gv2); // find the normal = (x, y, z) = (East, North, Up)
					buf_normals.put(normal);
					buf_normals.put(normal);
					buf_normals.put(normal);
					gt.setNormal(normal);
				}
				else if (usePerVertexNormal)
				{
					buf_normals.put(Geometry.computeVertexNormal(gv0, mapping));
					buf_normals.put(Geometry.computeVertexNormal(gv1, mapping));
					buf_normals.put(Geometry.computeVertexNormal(gv2, mapping));
				}
				
				triangleIndex++;
			}
		}
		buf_normals.flip();

        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Normals);
        glBufferData(GL_ARRAY_BUFFER, buf_normals, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
		System.out.println(" done!");
	}
	
	public void regenerate_normal_lines()
	{
		System.out.print("Regenerating normals lines...");
		
		final int numberTriangles = mapping.triangles.length;
		_renderToRender = false;
		
		// ensure that we de-reference these buffers
		buf_normal_lines = null;
		
		// create the buffers and maps
		buf_normal_lines = BufferUtils.createFloatBuffer(numberTriangles * VERTEX_SIZE * EDGE_SIZE); // 3 floats per 1 vertex, 2 vertices per normal line
		
		int triangleIndex = 0;
    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		while (iter.hasNext())
		{
			Triangle_dt dt = iter.next();
			Point_dt v0 = dt.p1();
			Point_dt v1 = dt.p2();
			Point_dt v2 = dt.p3();
			if (v0 != null && v1 != null && v2 != null)
			{ 
				// first get all the key representations of the vertices
				GeoTriangle gt = mapping.triangles[triangleIndex];

				// fill the vertex buffer with the correct vertices
				GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
				GeoVertex gv1 = mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
				GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
				
				float[] normal = Geometry.computeSurfaceNormal(gv0, gv1, gv2);
				
				// compute the normal lines
				float[] triangleCentre = new float[3];
				triangleCentre[0] = (float)((gv0.x() + gv1.x() + gv2.x())/3);
				triangleCentre[1] = (float)((gv0.y() + gv1.y() + gv2.y())/3);
				triangleCentre[2] = (float)((gv0.z() + gv1.z() + gv2.z())/3);
				buf_normal_lines.put(triangleCentre);
				float[] normal_line = Geometry.vectorPlus(triangleCentre, Geometry.vectorNormalize(normal)); 
				buf_normal_lines.put(normal_line);
				
				triangleIndex++;
			}
		}
		buf_normal_lines.flip();

        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_NormalLines);
        glBufferData(GL_ARRAY_BUFFER, buf_normal_lines, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
		System.out.println(" done!");
	}
	
	// depending on the mode, we update the color array
	public void repaint_triangulation()
	{
		System.out.print("Repainting mesh...");
	    switch(_mode)
		{
			case STRIKE_ONLY:
				repaint_StrikeNDip(true, false);
				break;
			case DIP_ONLY:
				repaint_StrikeNDip(false, true);
				break;
			case STRIKE_DIP:
				repaint_StrikeNDip(true, true);
				break;
				
			case STRIKE_DIP_SEARCH:
			case RANGE_QUERY:
			case PEAK_VALLEY_SADDLE:
			case ROUGHNESS:
				repaint_PurpleOverlay();
				break;
				
			case CURVATURE_GAUSSIAN:				
			case CURVATURE_MEAN:
			case CURVATURE_ANISOTROPIC:
			case CURVATURE_ANISOTROPIC_AVG:
			case ROUGHNESS_MAP:
			case CURVATURE_RING_MAPS:
				repaint_CurvatureMaps();
				break;
				
			case NONE:
			default:
				repaint_default();
				break;
		}
	    System.out.println(" done!");
	}
	
	public void repaint_default()
	{
		final int numberTriangles = mapping.triangles.length;
		_renderToRender = false;
		
		// ensure that we de-reference these buffers
		buf_colors = null;
		
		// recreate the color buffer
		buf_colors = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 color per 1 vertex
		
    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		while (iter.hasNext())
		{
			Triangle_dt dt = iter.next();
			Point_dt v0 = dt.p1();
			Point_dt v1 = dt.p2();
			Point_dt v2 = dt.p3();
			if (v0 != null && v1 != null && v2 != null) // we only care about the triangles with 3 points
			{
				buf_colors.put(CLRConst.GREEN);
				buf_colors.put(CLRConst.GREEN);
				buf_colors.put(CLRConst.GREEN);				
			}
		}
		buf_colors.flip();

        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Colors);
        glBufferData(GL_ARRAY_BUFFER, buf_colors, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        
		_renderToRender = true;
	}
	
	public void repaint_PurpleOverlay()
	{
		final int numberTriangles = mapping.triangles.length;
		_renderToRender = false;
		
		// ensure that we de-reference these buffers
		buf_colors = null;
		
		// recreate the color buffer
		buf_colors = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 color per 1 vertex
		
		int triangleIndex = 0;
    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		while (iter.hasNext())
		{
			Triangle_dt dt = iter.next();
			Point_dt v0 = dt.p1();
			Point_dt v1 = dt.p2();
			Point_dt v2 = dt.p3();
			if (v0 != null && v1 != null && v2 != null)
			{
				GeoTriangle gt = mapping.triangles[triangleIndex];
				
				if (gt.Selected)
				{
					buf_colors.put(CLRConst.PURPLE);
					buf_colors.put(CLRConst.PURPLE);
					buf_colors.put(CLRConst.PURPLE);					
				}
				else
				{
					buf_colors.put(CLRConst.GREEN);
					buf_colors.put(CLRConst.GREEN);
					buf_colors.put(CLRConst.GREEN);
				}
				
				triangleIndex++;
			}
		}
		buf_colors.flip();

        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Colors);
        glBufferData(GL_ARRAY_BUFFER, buf_colors, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        
		_renderToRender = true;
	}
	
	// regenerate the colors float buffer strike and dip colors
	public void repaint_StrikeNDip(boolean paintStrike, boolean paintDip)
	{
		final int numberTriangles = mapping.triangles.length;
		_renderToRender = false;
		
		// ensure that we de-reference these buffers
		buf_colors = null;
		
		// recreate the color buffer
		buf_colors = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 color per 1 vertex
		
		int triangleIndex = 0;
    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		while (iter.hasNext())
		{
			Triangle_dt dt = iter.next();
			Point_dt v0 = dt.p1();
			Point_dt v1 = dt.p2();
			Point_dt v2 = dt.p3();
			if (v0 != null && v1 != null && v2 != null)
			{
				GeoTriangle gt = mapping.triangles[triangleIndex];
				
				float[] clr = CLRConst.GREEN;
				
				if (perVertexStrikeDip)
				{
					for (Integer vk: gt.vertexkeys)
					{
						GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);
						
						double[] strike_n_dip = {gv.avgStrike, gv.avgDip};
						gt.StrikeDeg = strike_n_dip[0];
						gt.DipDeg = strike_n_dip[1];
						
						
						if (paintStrike && paintDip) clr = CLRConvert.GetStrikeDipColor(strike_n_dip[0], strike_n_dip[1]);
						else if (paintStrike) clr = CLRConvert.GetStrikeOnlyColor(strike_n_dip[0]);
						else if (paintDip) clr = CLRConvert.GetDipOnlyColor(strike_n_dip[1]);
						
						buf_colors.put(clr);
					}
				}
				else if (perTriangleStrikeDip)
				{
					GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
					GeoVertex gv1 = mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
					GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
					
					double[] strike_n_dip = Geometry.computeStrikeNDipAngle2(Geometry.computeSurfaceNormal(gv0, gv1, gv2), dLaserStrike, dLaserDip);
					gt.StrikeDeg = strike_n_dip[0];
					gt.DipDeg = strike_n_dip[1];
					
					if (paintStrike && paintDip) clr = CLRConvert.GetStrikeDipColor(strike_n_dip[0], strike_n_dip[1]);
					else if (paintStrike) clr = CLRConvert.GetStrikeOnlyColor(strike_n_dip[0]);
					else if (paintDip) clr = CLRConvert.GetDipOnlyColor(strike_n_dip[1]);
					
					buf_colors.put(clr);
					buf_colors.put(clr);
					buf_colors.put(clr);
				}

				triangleIndex++;
			}
		}
		buf_colors.flip();

        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Colors);
        glBufferData(GL_ARRAY_BUFFER, buf_colors, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        
		_renderToRender = true;
	}
	
	
	public void repaint_CurvatureMaps()
	{
		final int numberTriangles = mapping.triangles.length;
		_renderToRender = false;
		
		// ensure that we de-reference these buffers
		buf_colors = null;
		
		// recreate the color buffer
		buf_colors = BufferUtils.createFloatBuffer(numberTriangles * TRIANGLE_SIZE * VERTEX_SIZE); // 3 floats per 1 color per 1 vertex
		
		double[] colorRange = getColorScaleRange();
		double min = colorRange[0];
		double max = colorRange[1];
		
		int triangleIndex = 0;
    	Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		while (iter.hasNext())
		{
			Triangle_dt dt = iter.next();
			Point_dt v0 = dt.p1();
			Point_dt v1 = dt.p2();
			Point_dt v2 = dt.p3();
			if (v0 != null && v1 != null && v2 != null)
			{
				GeoTriangle gt = mapping.triangles[triangleIndex];

				GeoVertex[] vertices = {mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL),
										mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL),
										mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL) };

				for (GeoVertex gv: vertices)
				{
					double value = 0.0; // value is in the data range, we need to add the data_offset to it to place it within the color range
					double wave = 0.0;
					
					if (isMode_RoughnessMaps()) value = gv.curvatureRoughness;
					else if (isMode_CurvatureAnisotropic()) value = gv.curvatureAnisotropic;
					else if (isMode_CurvatureAnisotropicAvg()) value = gv.curvatureAnisotropicAveraged;
					else if (isMode_CurvatureMean()) value = gv.curvatureMean;
					else if (isMode_CurvatureGaussian()) value = gv.curvatureGaussian;
					else if (isMode_CurvatureRingMap()) value = gv.curvatureAvgRing;
					
					float[] rgb;
					if (isScaleMode_Log()) value = Math.log(value + data_offset);
					else if (isScaleMode_Sqrt()) value = Math.sqrt(value + data_offset);
					else value = value + data_offset;
					
					wave = CLRConvert.ConvertToWaveLength(min, max, value);
					rgb = CLRConvert.ConvertWaveLengthToRGB(wave);
					
					rgb[0] /= 255.0f;
					rgb[1] /= 255.0f;
					rgb[2] /= 255.0f;
					
					buf_colors.put(rgb);
				}

				triangleIndex++;
			}
		}

		buf_colors.flip();

        glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Colors);
        glBufferData(GL_ARRAY_BUFFER, buf_colors, GL_STATIC_DRAW);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        
		_renderToRender = true;
	}
	
	public void cleanup()
	{
		glDeleteBuffers(vboHandle_Vertices);
		glDeleteBuffers(vboHandle_Normals);
		glDeleteBuffers(vboHandle_Colors);

		_renderToRender = false;
	}
	
	private void render_mesh()
	{
		// render using VBO
		if (_renderToRender)
		{
	        glClear(GL_COLOR_BUFFER_BIT);

            glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Vertices);
            glVertexPointer(3, GL_FLOAT, 0, 0L); // vertex size = 3 => 3 floats per vertex
            glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Normals);
            glNormalPointer(GL_FLOAT, 0, 0L);
            glBindBuffer(GL_ARRAY_BUFFER, vboHandle_Colors);
            glColorPointer(3, GL_FLOAT, 0, 0L); // color size = 3 => 3 floats per color
            
            glEnableClientState(GL_VERTEX_ARRAY);
            glEnableClientState(GL_NORMAL_ARRAY);
            glEnableClientState(GL_COLOR_ARRAY);
            
            glDrawArrays(GL_TRIANGLES, 0, mapping.triangles.length*3);
            
            glDisableClientState(GL_COLOR_ARRAY);
            glDisableClientState(GL_NORMAL_ARRAY);
            glDisableClientState(GL_VERTEX_ARRAY);
            
            
            // display the vector normals
            if (showNormals)
            {            	
                glBindBuffer(GL_ARRAY_BUFFER, vboHandle_NormalLines);
                glVertexPointer(3, GL_FLOAT, 0, 0L); // vertex size = 3 => 3 floats per vertex
                
            	glColor3f(1.0f, 1.0f, 1.0f);
            	glEnableClientState(GL_VERTEX_ARRAY);
            	glDrawArrays(GL_LINES, 0, mapping.triangles.length*2);
            	glDisableClientState(GL_VERTEX_ARRAY);
            }
            
            
//            glBindBuffer(GL_ARRAY_BUFFER, 0); // TODO: find out if this line is truly un-needed
		}	
	}
	
	public void render()
	{
		if (_wireframe) glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
		else if (_point) glPolygonMode(GL_FRONT_AND_BACK, GL_POINT);
		else glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
		
		// ensure we are ready to render before we alter the rendering data!
		if (_renderToRender)
		{
			if (updateComplete)
			{
				regenerate_complete();
				updateComplete = false;
			}
			
			if (updateGeometry)
			{
				regenerate_geometry();
				updateGeometry = false;
			}
			
			if (updateColors)
			{
				repaint_triangulation();
				updateColors = false;
			}
			
			if (usePerTriangleNormal)
			{
				regenerate_normals();
				usePerTriangleNormal = false;
			}
			
			if (usePerVertexNormal)
			{
				regenerate_normals();
				usePerVertexNormal = false;
			}
		}

		glPushMatrix();
		
		float x = (float)(xBoundMin + getBoundX_Legnth()/2);
		float y = (float)(yBoundMin + getBoundY_Legnth()/2);
		float z = (float)(zBoundMin);
		
		glScalef(MESH_ZOOM, MESH_ZOOM, 1.0f);
		glTranslated(x, y, z);
		glRotatef(_zRot, 0.0f, 0.0f, 1.0f);
		glRotatef(_yRot, 0.0f, 1.0f, 0.0f);
		glRotatef(_xRot, 1.0f, 0.0f, 0.0f);
		glTranslated(-x, -y, -z);
		
		render_mesh();
		
		glPopMatrix();
	}
	
	// TODO: this is legacy code. only computeRoughnessV1 uses this method
	// computeRoughnessV1 requires that each triangle know it's neighbouring triangles thus we should add the indices to the DelaunyMapping
	// inorder to remove this
	private float[] computeSurfaceNormal(Point_dt v0, Point_dt v1, Point_dt v2)
	{
		float[] normal = {0.0f, 0.0f, 0.0f};
		Vector p0 = new Vector(new double[] {v0.x(), v0.y(), v0.z()});
		Vector p1 = new Vector(new double[] {v1.x(), v1.y(), v1.z()});
		Vector p2 = new Vector(new double[] {v2.x(), v2.y(), v2.z()});
		p0.enableImmutable(); p1.enableImmutable(); p2.enableImmutable();
		
		Vector U = p1.minus(p0);
		Vector V = p2.minus(p0);
		normal[0] = (float)((U.get(1) * V.get(2)) - (U.get(2) * V.get(1)));
		normal[1] = (float)((U.get(2) * V.get(0)) - (U.get(0) * V.get(2)));
		normal[2] = (float)((U.get(0) * V.get(1)) - (U.get(1) * V.get(0)));
		
		return normal;
	}
	
	// computes the strike && dip angle given a normal = (x, y, z) = (East, North, Up) of a plane
	// returns (strike, dip)
	private double[] computeStrikeNDipAngle(float[] n)
	{
		double[] angles = new double[2];
		
//		// *** in theory our normals SHOULD all point in the right direction => in the direction of the camera head ***
//		// flip the sign of the normal vector if the normal points down
		if (n[2] < 0)
		{
			n[0] = -n[0];
			n[1] = -n[1];
			n[2] = -n[2];
		}
		
		// strike = -tan^-1(y/x)
		// dip = tan^-1  (sqrt(x^2 + y^2)/z)
		
//		double strike = Math.atan(n[1]/n[0]);
//		double strikeDeg = Math.toDegrees(strike);
//		
//		double dip = Math.atan((Math.sqrt(n[1]*n[1] + n[0]*n[0])/n[2]));
//		double dipDeg = Math.toDegrees(dip);
		
		
//		// orignal strike n  dip code
//		double strike = -Math.atan(n[0]/n[1]);
		double strike = Math.atan(n[1]/n[0]); // this is what it should be (North/East) but it does not produce a good result....
		double strikeDeg = Math.toDegrees(strike);
		
		double dip = Math.atan((Math.sqrt(n[0]*n[0] + n[1]*n[1])/n[2]));
		double dipDeg = Math.toDegrees(dip);
		
		// go round and round until we are in range
		strikeDeg = strikeDeg - 90.0;
		while (strikeDeg >= 360) strikeDeg = strikeDeg - 360.0;
		while (strikeDeg < 0) strikeDeg = strikeDeg + 360.0;
		
		angles[0] = strikeDeg;
		angles[1] = dipDeg;
		return angles;
	}
	
	// find the angles between the normals of the incidient triangles and sum up the total angle
	private void computeRoughness()
	{
		int triangleIndex = 0;
		Iterator<Triangle_dt> iter = _delaunay.trianglesIterator();
		
		while(iter.hasNext())
		{
			Triangle_dt t = iter.next();
			if (t.p1() != null && t.p2() != null && t.p3() != null) // we only want triangles with at least 3 vertices
			{
				Triangle_dt[] neighbours = {t.next_12(), t.next_23(), t.next_31()};
				float[] t_norm = computeSurfaceNormal(t.p1(), t.p2(), t.p3());
				double t_magnitude = t_norm[0]*t_norm[0] + t_norm[1]*t_norm[1] + t_norm[2]*t_norm[2];

				double roughness = 0.0;
				for (Triangle_dt dt: neighbours)
				{
					Point_dt v0 = dt.p1();
					Point_dt v1 = dt.p2();
					Point_dt v2 = dt.p3();
					
					// find the angles between the normals
					if (v0 != null && v1 != null && v2 != null)
					{
						float[] dt_norm = computeSurfaceNormal(v0, v1, v2);
						double dt_magnitude = dt_norm[0]*dt_norm[0] + dt_norm[1]*dt_norm[1] + dt_norm[2]*dt_norm[2];
						double dot = t_norm[0]*dt_norm[0] + t_norm[1]*dt_norm[1] + t_norm[2]*dt_norm[2];
						double tmp = Math.acos(dot/(t_magnitude*dt_magnitude));
						if (!Double.isNaN(tmp)) roughness += tmp; 
					}
				}
				mapping.triangles[triangleIndex].setRoughness(roughness);
				triangleIndex++;
			}
		}
	}
	
	// label all the vertices as a peak, valley or saddle
	public void computePeakValleySaddle()
	{
		for (int i = 0; i < mapping.vertices.length; i++)
		{
			int peaks = 0;
			GeoVertex gv = mapping.vertices[i];
			for (Integer vKey: gv.neighbourkeys)
			{
				GeoVertex n = mapping.vertices[vKey];
				if (n.z() < gv.z()) peaks++;
				else if (n.z() == gv.z())
				{
					if (n.y() < gv.y()) peaks++;
					else if (n.y() == gv.y())
					{
						if (n.x() < gv.x()) peaks++;
					}
				}
			}
			gv.setPeaks(peaks);
		}
	}
	
	// highlight all triangles within some range
	public void searchRange(double x, double y, double w, double h)
	{
		if (_tree == null) // 
		{
			ArrayList<PointND> found = _tree.range_find(new QuadBoundingBox(y, x, y+h, x+w), _tree.getRoot());

//			// reset all the selected triangles
//			_numOfSelectedTriangles = 0;
//			for(GeoTriangle gi: _triangleMap.values()) gi.Selected = false;
//
//			for (PointND p: found)
//			{
//				String key = p.get_coord(0)+","+p.get_coord(1)+","+p.get_coord(2);
//				if (_vertexMap.containsKey(key))
//				{
//					HashSet<String> triangleKeys = _vertexMap.get(key).getTriangleKeys();
//					for (String tk: triangleKeys)
//					{
//						if (_triangleMap.containsKey(tk))
//						{
//							_triangleMap.get(tk).Selected = true;
//							_numOfSelectedTriangles++;
//						}
//					}
//				}
//			}
			
			_numOfSelectedTriangles = 0;
			
			for (GeoTriangle gt: mapping.triangles) gt.Selected = false;
			
			
			System.out.println("Search Range Query ==> Found " + found.size() + " vertices within range ("+x+","+y+") -- ("+(x+w)+","+(y+h)+")");
		}
		else
		{
			System.out.println("Quad-Tree is NULL, aborting range query...");
		}
	}
	
	// search the triangulation using mouse dragging
	public void searchRange_dragWindow(float x1, float y1, float x2, float y2)
	{
		// TODO: do proper projection before doing a range query
		float x = x1*UIConst.POINT_CLOUD_SCALING;
		float y = y1*UIConst.POINT_CLOUD_SCALING;
		float w = Math.abs(x1 - x2)*UIConst.POINT_CLOUD_SCALING;
		float h = Math.abs(y1 - y2)*UIConst.POINT_CLOUD_SCALING;
		searchRange(x, y, w, h);
	}
	
	// compute the strike and dip per 1-ring vertex and per triangle
	public void computeStrikeDip_All()
	{
		// average the strike and dips per vertex?
		System.out.print("Computing strike n dip per vertex... ");
		for (GeoVertex v: mapping.vertices)
		{
			GeoVertex gv = mapping.getVertex(v.id, CURRENT_SMOOTHNESS_LEVEL);
			HashSet<Integer> visited = new HashSet<Integer>();
			ArrayList<ArrayList<GeoVertex>> rings = findRingNeighbourhood(gv, 1, visited);
			double avgStrike = 0.0;
			double avgDip = 0.0;
			
			int num = 0;
			for (ArrayList<GeoVertex> r: rings)
			{
				for (GeoVertex found: r)
				{
					double[] snd = Geometry.computeStrikeNDipAngle2(Geometry.computeVertexNormal(found, mapping), dLaserStrike, dLaserDip);
					avgStrike += snd[0];
					avgDip += snd[1];
					num++;
				}
			}
			
			avgStrike /= num;
			avgDip /= num;
			gv.avgStrike = avgStrike;
			gv.avgDip = avgDip;
		}
		System.out.println(" done!");
		
		System.out.print("Computing strike n dip per triangle... ");
		for (int tk = 0; tk < mapping.triangles.length; tk++)
		{
			GeoTriangle gt = mapping.triangles[tk];
			GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv1 = mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
			double[] strike_dip = Geometry.computeStrikeNDipAngle2(Geometry.computeSurfaceNormal(gv0, gv1, gv2), dLaserStrike, dLaserDip);
			gt.StrikeDeg = strike_dip[0];
			gt.DipDeg = strike_dip[1];
		}
		System.out.println(" done!");
	}
	
	public void searchStrikeRange(double strike_from, double strike_to)
	{
		_numOfSelectedTriangles = 0;
		for (int tk = 0; tk < mapping.triangles.length; tk++)
		{
			GeoTriangle gt = mapping.triangles[tk];
			gt.Selected = false;
			
			GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv1= mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
			double[] strikeDip = Geometry.computeStrikeNDipAngle2(Geometry.computeSurfaceNormal(gv0, gv1, gv2), dLaserStrike, dLaserDip);
			
			if (strikeDip[0] >= strike_from && strikeDip[0] <= strike_to) // we only care about the strike for now
			{
				gt.Selected = true;
				_numOfSelectedTriangles++;
			}
		}
		updateColors = true;
		System.out.println("Search Strike ==> Found " + _numOfSelectedTriangles + " matches in " + mapping.triangles.length + " triangles.");
	}
	
	public void searchDipRange(double dip_from, double dip_to)
	{
		_numOfSelectedTriangles = 0;
		for (int tk = 0; tk < mapping.triangles.length; tk++)
		{
			GeoTriangle gt = mapping.triangles[tk];
			gt.Selected = false;
			
			GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv1= mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
			double[] strikeDip = Geometry.computeStrikeNDipAngle2(Geometry.computeSurfaceNormal(gv0, gv1, gv2), dLaserStrike, dLaserDip);
			
			if (strikeDip[1] >= dip_from && strikeDip[1] <= dip_to) // we only care about the dip for now

			{
				gt.Selected = true;
				_numOfSelectedTriangles++;
			}
		}
		updateColors = true;
		System.out.println("Search Dip ==> Found " + _numOfSelectedTriangles + " matches in " + mapping.triangles.length + " triangles.");
	}
	
	public void searchStrikeDipRange(double strike_from, double strike_to, double dip_from, double dip_to)
	{
		_numOfSelectedTriangles = 0;
		for (int tk = 0; tk < mapping.triangles.length; tk++)
		{
			GeoTriangle gt = mapping.triangles[tk];
			gt.Selected = false;
			
			GeoVertex gv0 = mapping.getVertex(gt.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv1= mapping.getVertex(gt.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex gv2 = mapping.getVertex(gt.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
			double[] strikeDip = Geometry.computeStrikeNDipAngle2(Geometry.computeSurfaceNormal(gv0, gv1, gv2), dLaserStrike, dLaserDip);
			
			if (strikeDip[0] >= strike_from && strikeDip[0] <= strike_to && strikeDip[1] >= dip_from && strikeDip[1] <= dip_to) // we only care about the dip for now
			{
				gt.Selected = true;
				_numOfSelectedTriangles++;
			}
		}
		updateColors = true;
		System.out.println("Search Strike & Dip ==> Found " + _numOfSelectedTriangles + " matches in " + mapping.triangles.length + " triangles.");
	}
	
	
	// find all triangles that has some "roughness"
	public void searchRoughness(double rough_from, double rough_to)
	{
		double MAX_ROUGHNESS = 540;
		
		_numOfSelectedTriangles = 0;
		for (int tk = 0; tk < mapping.triangles.length; tk++)
		{
			GeoTriangle gv = mapping.triangles[tk];
			double percentage = (gv.RoughnessDeg/MAX_ROUGHNESS) * 100;
			boolean bInRange = percentage >= rough_from && percentage <= rough_to;
			gv.Selected = false;
			if (bInRange)
			{
				gv.Selected = true;
				_numOfSelectedTriangles++;
			}
		}
		System.out.println("Search Roughness v1 ==> Found " + _numOfSelectedTriangles + " matches in " + mapping.triangles.length + " triangles.");
	}
	
	public void searchRoughnessRatio(double from, double to)
	{
		_numOfSelectedTriangles = 0;
		for (int tk = 0; tk < mapping.triangles.length; tk++) mapping.triangles[tk].Selected = false; // reset the selection flag
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);
			if (gv.getPeakRatio() >= from && gv.getPeakRatio() <= to)
			{
				for (Integer key: gv.trianglekeys)
				{
					GeoTriangle gt = mapping.triangles[key];
					if (!gt.Selected)
					{
						gt.Selected = true;
						_numOfSelectedTriangles++;
					}					
				}
			}
		}
		updateColors = true;
		System.out.println("Search Peak Ratio ==> Found " + _numOfSelectedTriangles + " matches in " + mapping.triangles.length + " triangles.");
	}
	
	public ArrayList<GeoEdge> findEdgeSubset(GeoVertex src, GeoVertex cur, Integer curKey, HashSet<Integer> visited, double radius)
	{
		ArrayList<GeoEdge> ret = new ArrayList<GeoEdge>();
		
		// recursive stopping condition
		// can't add the neighbour keys like in findVertexSubset as we'll skip edges due to having visited the vertex
		// but not having checked to see if the edge has been visited by recursion
		if (visited.contains(curKey)) return ret;
		visited.add(curKey);
		
		for (Integer nk: cur.neighbourkeys)
		{
			GeoVertex node = mapping.getVertex(nk, CURRENT_SMOOTHNESS_LEVEL);
			GeoEdge e = new GeoEdge(cur, node);
			
			double dist = src.distance3D(node);
			if (dist < radius)
			{
				ret.add(e);
				ret.addAll(findEdgeSubset(src, node, nk, visited, radius));
			}
			else
			{
				double diff = dist - radius;
				if (diff == 0.0 || diff < 0.0001)
				{
					ret.add(e);
				}
				else
				{
					double param = Geometry.SphereLineIntersection(src, e.start, e.end, radius);
					GeoVertex gv = e.parametricGeoVertex(param);  // TODO: abstract the parametricGeoVertex method or create the map of edges
					GeoEdge e_boundary = new GeoEdge(e.start, gv, e.end); // edge that should go right up to the circle radius, special parametic edge
					ret.add(e_boundary);
				}
			}
		}
		return ret;
	}
	
	public ArrayList<GeoVertex> findVertexSubset(GeoVertex src, GeoVertex cur, double radius, HashSet<Integer> visited)
	{
		ArrayList<GeoVertex> ret = new ArrayList<GeoVertex>();
		for (Integer nk: cur.neighbourkeys)
		{
			if (visited.contains(nk)) continue;
			visited.add(nk);

			GeoVertex node = mapping.getVertex(nk, CURRENT_SMOOTHNESS_LEVEL);
			double dist = src.distance3D(node);
			if (dist < radius)
			{
				ret.add(node);
//				ret.addAll(findVertexSubset(src, mapping.vertices[nk], radius, visited));
				ret.addAll(findVertexSubset(src, node, radius, visited));
			}
			else
			{
				GeoEdge e = new GeoEdge(cur, node);
				GeoVertex gv = e.parametricGeoVertex(radius/dist);  // TODO: abstract the parametricGeoVertex method or create the map of edges
				ret.add(gv);
			}
		}
		
		return ret;
	}

	public ArrayList<GeoVertex> findEdgeVertexSubset(GeoVertex src, GeoVertex cur, double radius, HashSet<Integer> visited)
	{
		ArrayList<GeoVertex> ret = new ArrayList<GeoVertex>();
		for (Integer nk: cur.neighbourkeys)
		{
			if (visited.contains(nk)) continue;
			visited.add(nk);

			GeoVertex node = mapping.getVertex(nk, CURRENT_SMOOTHNESS_LEVEL);
			double dist = src.distance3D(node);
			if (dist <= radius)
			{
				if (dist == radius) ret.add(node);
//				else ret.addAll(findEdgeVertexSubset(src, mapping.vertices[nk], radius, visited));
				else ret.addAll(findEdgeVertexSubset(src, node, radius, visited));
			}
			else
			{
				GeoEdge e = new GeoEdge(cur, node);
				double param = Geometry.SphereLineIntersection(src, e.start, e.end, radius);
				GeoVertex gv = e.parametricGeoVertex(param);
				ret.add(gv);
			}
		}
		
		return ret;
	}
	
	// find all the edges that exist within a local window defined around src with radius r
	public ArrayList<GeoEdge> localWindow_Edges(GeoVertex src, double r)
	{
		ArrayList<GeoEdge> edges = new ArrayList<GeoEdge>();
		HashSet<Integer> visited = new HashSet<Integer>();
		Integer srcIndex = src.id;
		edges = findEdgeSubset(src, src, srcIndex, visited, r);
		return edges;
	}
	
	// find all the vertices that exist within a local window defined around src with radius r
	public ArrayList<GeoVertex> localWindow_Vertices(GeoVertex src, double r)
	{
		ArrayList<GeoVertex> vertices = new ArrayList<GeoVertex>();
		HashSet<Integer> visited = new HashSet<Integer>();
		vertices = findVertexSubset(src, src, r, visited);
		return vertices;
	}
	
	// find all the EDGE vertices that exist within a local window defined around src with radius r
	public ArrayList<GeoVertex> localWindow_EdgeVertices(GeoVertex src, double r)
	{
		ArrayList<GeoVertex> vertices = new ArrayList<GeoVertex>();
		HashSet<Integer> visited = new HashSet<Integer>();
		vertices = findEdgeVertexSubset(src, src, r, visited);
		return vertices;		
	}

	// find the n-ring neighbour hood centered on the src vertex
	public ArrayList<ArrayList<GeoVertex>> findRingNeighbourhood(GeoVertex src, int n, HashSet<Integer> visited)
	{
		// create the list which will contain all the rings
		ArrayList<ArrayList<GeoVertex>> rings = new ArrayList<ArrayList<GeoVertex>>();
		for (int i = 0; i < n+1; i++) rings.add(new ArrayList<GeoVertex>());
		
		LinkedList<GeoVertex> current = new LinkedList<GeoVertex>(); // current i-ring
		LinkedList<GeoVertex> next = new LinkedList<GeoVertex>(); // possible vertices for (i+1)-ring
		
		visited.add(src.id); // init the queue
		current.add(src);
		rings.get(0).add(src);

		for (int i = 1; i < n+1; i++)
		{
			while (current.isEmpty() == false)
			{
				GeoVertex g = mapping.getVertex(current.removeFirst().id, CURRENT_SMOOTHNESS_LEVEL);
				for (Integer nk: g.neighbourkeys)
				{
					if (visited.contains(nk)) continue; // skip all vertices that we have visited at any ring level
					visited.add(nk);

					GeoVertex gv = mapping.getVertex(nk, CURRENT_SMOOTHNESS_LEVEL);
					rings.get(i).add(gv);
					next.add(gv);
				}
			}
			current = next;
			next = new LinkedList<GeoVertex>();
		}
		return rings;
	}
	
	// reset each vertices smoothed points list
	public void smoothMesh_Reset()
	{
		for (GeoVertex g: mapping.vertices) g.clearSmoothedVertices();
		CURRENT_SMOOTHNESS_LEVEL = 0;
	}
	
	public void smoothMesh_geometricLaplacian(int iterations)
	{	
		smoothMesh_Reset();
		MAX_SMOOTHNESS_LEVEL = iterations;
		System.out.println("Performing " + iterations +" iterations of Geometric Laplacian smoothing... ");
		for (int i = 0; i < iterations; i++)
		{
			boolean bUseSmoothedVertex = i > 0 || CURRENT_SMOOTHNESS_LEVEL > 0;
			for (GeoVertex g: mapping.vertices)
			{
				GeoVertex smoothed = null;
				if (bUseSmoothedVertex) smoothed = Smoothing.geometricLaplacian(g.getSmoothestVertex(), g.neighbourkeys, i, mapping);
				else smoothed = Smoothing.geometricLaplacian(g, g.neighbourkeys, i, mapping);
				g.addLOSVertex(smoothed);
			}
			System.out.println("... completed iteration " + (i+1) + " ...");
		}
	}
	
	public void smoothMesh_Adaptive(int iterations, double scale)
	{
		smoothMesh_Reset();
		MAX_SMOOTHNESS_LEVEL = iterations;
		double radius = 2.0*scale*(xBoundMax - xBoundMin);
		
		System.out.println("Performing " + iterations +" iterations of Adaptive smoothing with radius="+radius+"... ");
		for (int i = 0; i < iterations; i++)
		{
			boolean bUseSmoothedVertex = i > 0 || CURRENT_SMOOTHNESS_LEVEL > 0;
			for (GeoVertex g: mapping.vertices)
			{
				GeoVertex smoothed  = null;
				GeoVertex gv = null;
				
				if (bUseSmoothedVertex) gv = g.getSmoothestVertex();
				else gv = g;
				
				ArrayList<GeoVertex> edge_points = localWindow_EdgeVertices(gv, radius);
				smoothed = Smoothing.adaptiveSmoothing(gv, edge_points);
				g.addLOSVertex(smoothed);
			}
			System.out.println("... completed iteration " + (i+1) + " ...");
		}
	}
	
	public void smoothMesh_Taubin(int iterations)
	{
		smoothMesh_Reset();
		MAX_SMOOTHNESS_LEVEL = iterations;
		double lambda = 0.6307;
		double mu = -0.6732;
		
		System.out.println("Performing " + iterations +" iterations of Taubin smoothing ... ");
		for (int i = 0; i < iterations; i++)
		{
			boolean bUseSmoothedVertex = i > 0 || CURRENT_SMOOTHNESS_LEVEL > 0;
			for (GeoVertex g: mapping.vertices)
			{
				Point_dt laplacian = null;
				Point_dt current = null;
				if (bUseSmoothedVertex)
				{
					laplacian = Smoothing.geometricLaplacian(g.getSmoothestVertex(), g.neighbourkeys, i, mapping);
					current = g.getSmoothestVertex();
				}
				else
				{
					laplacian = Smoothing.geometricLaplacian(g, g.neighbourkeys, i, mapping);
					current = g;
					
				}

				double[] d = new double[3];
				d[0] = Math.pow((current.x() - lambda*laplacian.x())*(current.x() - mu*laplacian.x()), i / 2);
				d[1] = Math.pow((current.y() - lambda*laplacian.x())*(current.y() - mu*laplacian.y()), i / 2);
				d[2] = Math.pow((current.z() - lambda*laplacian.x())*(current.z() - mu*laplacian.z()), i / 2);

				GeoVertex smoothed = new GeoVertex(new Point_dt(d[0], d[1], d[2]));
				smoothed.neighbourkeys = g.neighbourkeys;
				smoothed.trianglekeys = g.trianglekeys;
				smoothed.id = g.id;
				g.addLOSVertex(smoothed);
			}
			System.out.println("... completed iteration " + (i+1) + " ...");
		}
	}
	
	public void computeRoughnessMaps(double scale, int smoothingIterations)
	{
		System.out.println("Computing Roughness Map... ");
		MAX_SMOOTHNESS_LEVEL = smoothingIterations;
		
		// first smooth the mesh
		System.out.println("\tSmoothing original mesh... ");
		smoothMesh_Adaptive(smoothingIterations, scale);
		System.out.println("\tdone!");
		
		// save this to restore later, setting this to zero allows us to use the original GeoVertices
		CURRENT_SMOOTHNESS_LEVEL = MAX_SMOOTHNESS_LEVEL;

		// find the curvature values for the smoothed mesh
		System.out.println("\tComputing kmax curvature for smoothed mesh... ");
		computeAnisotropicCurvatureMaps(scale);
		System.out.println("\tdone!");
		
		// find the curvature values for the original mesh
		System.out.println("\tComputing kmax curvature for original mesh... ");
		CURRENT_SMOOTHNESS_LEVEL = 0;
		computeAnisotropicCurvatureMaps(scale);
		CURRENT_SMOOTHNESS_LEVEL = MAX_SMOOTHNESS_LEVEL;
		System.out.println("\tdone!");
		
		// now compute the average curvature for each vertex as the mean of the local window
		System.out.println("\tComputing kav curvature for smoothed mesh... ");
		computeAvgAnisotropicCurvatureMaps(2*scale); // double radius size allows us to detect noise that has bee surpressed
		System.out.println("\tdone!");
		
		System.out.println("\tComputing kav curvature for original mesh... ");
		CURRENT_SMOOTHNESS_LEVEL = 0;
		computeAvgAnisotropicCurvatureMaps(2*scale);
		CURRENT_SMOOTHNESS_LEVEL = MAX_SMOOTHNESS_LEVEL;
		System.out.println("\tdone!");
		
		double maxRough = Double.MIN_VALUE;
		double minRough = Double.MAX_VALUE;
		double avgRough = 0.0;
		
		// compute the asymmetric difference
		System.out.print("\tComputing asymmetric difference... ");
		for (GeoVertex gv: mapping.vertices)
		{
			double kav = gv.curvatureAnisotropicAveraged; // kav of original vertex
			double kav_s = gv.getSmoothestVertex().curvatureAnisotropicAveraged; // kav of smoothed vertex
			double roughness = kav > kav_s ? kav - kav_s : 0.0;
			gv.curvatureRoughness = roughness;
			
			maxRough = Math.max(roughness, maxRough);
			minRough = Math.min(roughness, minRough);
			avgRough += roughness;
		}
		setColorScaleRange(minRough, maxRough);
		System.out.println(" done!");
		
		avgRough /= mapping.vertices.length;
		System.out.println("Roughness[Max:"+maxRough+"  Min:"+minRough+"  Avg:"+avgRough+"]");
		
		CURRENT_SMOOTHNESS_LEVEL = 0;
		updateColors = true;
	}
	
	public void computeRoughnessMapV2(double scale, int smoothingIterations)
	{
		System.out.println("Computing Roughness Map... ");

		// first smooth the mesh
		System.out.println("\tSmoothing original mesh... ");
		smoothMesh_Adaptive(smoothingIterations, scale);
		System.out.println("\tdone!");
		
		// save this to restore later, setting this to zero allows us to use the original GeoVertices
		int actual_smoothness_lvl = smoothingIterations;
		CURRENT_SMOOTHNESS_LEVEL = 0; // compute principle curvatures for the non-smoothed mesh
		computePrincipleCurvatures_Taubin();
		
		CURRENT_SMOOTHNESS_LEVEL = actual_smoothness_lvl; // now for the smoothed mesh
		computePrincipleCurvatures_Taubin();
		
		CURRENT_SMOOTHNESS_LEVEL = 0; // now compute the average mean curvature for the non-smoothed mesh
		computeAvgMeanCurvatureMaps(2*scale);
		
		CURRENT_SMOOTHNESS_LEVEL = actual_smoothness_lvl; // again for the smoothed mesh
		computeAvgMeanCurvatureMaps(2*scale);
		
		
		double maxRough = Double.MIN_VALUE;
		double minRough = Double.MAX_VALUE;
		
		// compute the asymmetric difference
		System.out.print("\tComputing asymmetric difference... ");
		for (GeoVertex gv: mapping.vertices)
		{
			double kav = gv.curvatureMean; // kav of original vertex
			double kav_s = gv.getSmoothestVertex().curvatureMeanAveraged; // kav of smoothed vertex
			double roughness = kav > kav_s ? kav - kav_s : 0.0;
			gv.curvatureRoughness = roughness;
			System.out.println(kav + " " + kav_s);
			
			maxRough = Math.max(roughness, maxRough);
			minRough = Math.min(roughness, minRough);
		}
		setColorScaleRange(minRough, maxRough);
		System.out.println(" done!");
		
		updateColors = true;
	}
	
	public void computeAvgMeanCurvatureMaps(double scale)
	{
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);
			double kav = 0.0;
			
			ArrayList<GeoVertex> window = localWindow_Vertices(gv, scale);
			for (GeoVertex g: window)
			{
				kav += g.curvatureMean;
			}
			
			kav /= window.size();
			
			gv.curvatureMeanAveraged = kav;
		}
	}
	
	public void computeAvgAnisotropicCurvatureMaps(double scale)
	{
		double max = Double.MIN_VALUE;
		double min = Double.MAX_VALUE;
		double avg = 0.0;
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);
			double kav = 0.0;
			
			ArrayList<GeoVertex> window = localWindow_Vertices(gv, scale);
			for (GeoVertex g: window)
			{
				kav += g.curvatureAnisotropic;
			}
			
			kav /= window.size();
			
			max = Math.max(kav, max);
			min = Math.min(kav, min);
			avg += kav;
			
			gv.curvatureAnisotropicAveraged = kav;
		}
		
		System.out.println("Averaged Anisotropic= max("+max+")  min("+min+")  avg("+(avg/mapping.vertices.length)+")");
	}
	
	public void computeAnisotropicCurvatureMaps(double scale)
	{
		double anisotropic_max_curvature = Double.MIN_VALUE;
		double anisotropic_min_curvature = Double.MAX_VALUE;
		double anisotropic_avg_curvature = 0.0;
		
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);
			
			// for each vertex we find a subset of edges which span the subset of vertices contained within the local window
			ArrayList<GeoEdge> subset = localWindow_Edges(gv, scale*(xBoundMax - xBoundMin));

			double[] anisotropic_curvature = Curvatures.computeAnisotropicCurvature(subset, mapping, CURRENT_SMOOTHNESS_LEVEL);
			
			// find kmax as it is the curvature value which represents roughness
			double kmax = Math.max(anisotropic_curvature[0], anisotropic_curvature[1]);
			kmax = Math.max(kmax, anisotropic_curvature[2]);
			
			if (anisotropic_max_curvature < kmax) anisotropic_max_curvature = kmax;
			if (anisotropic_min_curvature > kmax) anisotropic_min_curvature = kmax;
			if (Double.isNaN(kmax) == false)
			{
				anisotropic_avg_curvature += kmax;
				gv.curvatureAnisotropic = kmax;
			}
		}
		
		System.out.println("Anisotropic= max("+anisotropic_max_curvature+")  min("+anisotropic_min_curvature+")  avg("+(anisotropic_avg_curvature/mapping.vertices.length)+")");
	}

	public void computeMeanCurvatureMaps(double scale)
	{
		double mean_max_curvature = Double.MIN_VALUE;
		double mean_min_curvature = Double.MAX_VALUE;
		double mean_avg_curvature = 0.0;
		
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);
			
			// for each vertex we find a subset of edges which span the subset of vertices contained within the local window
			ArrayList<GeoEdge> subset = localWindow_Edges(gv, scale*(xBoundMax - xBoundMin));

			double mean_curvature = Curvatures.computeMeanCurvature(subset, mapping, CURRENT_SMOOTHNESS_LEVEL);
			if (mean_max_curvature < mean_curvature) mean_max_curvature = mean_curvature;
			if (mean_min_curvature > mean_curvature) mean_min_curvature = mean_curvature;
			if (Double.isNaN(mean_curvature) == false)
			{
				mean_avg_curvature += mean_curvature;
				gv.curvatureMean = mean_curvature;
			}
			else
			{
				gv.curvatureMean = 0.0;
			}
		}
		System.out.println("Mean= max("+mean_max_curvature+")  min("+mean_min_curvature+")  avg("+(mean_avg_curvature/mapping.vertices.length)+")");
		setColorScaleRange(mean_min_curvature, mean_max_curvature);
	}
	
	public void computeGaussianCurvatureMaps()
	{
		double gaussian_max_curvature = Double.MIN_VALUE;
		double gaussian_min_curvature = Double.MAX_VALUE;
		double gaussian_avg_curvature = 0.0;
		
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);

			double gaussian_curvature = Geometry.angleOfDefect(gv, mapping);
			if (gaussian_max_curvature < gaussian_curvature) gaussian_max_curvature = gaussian_curvature;
			if (gaussian_min_curvature > gaussian_curvature) gaussian_min_curvature = gaussian_curvature;
			if (Double.isNaN(gaussian_curvature) == false)
			{
				gaussian_avg_curvature += gaussian_curvature;
				gv.curvatureGaussian = gaussian_curvature;
			}
			else
			{
				gv.curvatureGaussian = 0.0;
			}
		}
		System.out.println("Gaussian= max("+gaussian_max_curvature+")  min("+gaussian_min_curvature+")  avg("+(gaussian_avg_curvature/mapping.vertices.length)+")");
		setColorScaleRange(gaussian_min_curvature, gaussian_max_curvature); 
	}
	
	// compute k1 and k2 via Tauins estimation method
	public void computePrincipleCurvatures_Taubin()
	{
		double pc_gaus_min = Double.MAX_VALUE;
		double pc_gaus_max = Double.MIN_VALUE;
		
		double pc_mean_min = Double.MAX_VALUE;
		double pc_mean_max = Double.MIN_VALUE;
		
		System.out.print("Computing Principle Curvatures (Taubin)... ");
		// first compute the vertex normal, shared triangle area (per neighbour per vertex) and total triangle area per vertex
		HashMap<Integer, HashMap<Integer, Double>> sharedTriangleAreas = new HashMap<Integer, HashMap<Integer,Double>>();
		HashMap<Integer, Double> totalTriangleAreas = new HashMap<Integer, Double>();
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL);
			float[] n = Geometry.computeVertexNormal(gv, mapping); // compute the normal for the vertex that will be rendered to the screen
			
			HashMap<Integer, Double> shared = new HashMap<Integer, Double>();
			double SWi = 0.0; // total area of the incident triangles around vertex v_i, we will need this value to scale everything to 1
			for (Integer nk: gv.neighbourkeys)
			{
				Integer tk_1 = null; // triangle keys that are in common with the neighbour
				Integer tk_2 = null;
				
				GeoVertex gvn = mapping.getVertex(nk, CURRENT_SMOOTHNESS_LEVEL);
				for (Integer tk: gvn.trianglekeys)
				{
					if (gv.trianglekeys.contains(tk))
					{
						if (tk_1 == null) tk_1 = tk;
						else if (tk_2 == null)
						{
							tk_2 = tk;
							break;
						}				
					}
				}
				
				double sharedArea = 0.0; // this is the w_ij value
				// add the triangle areas to the total surface area SWi
				if (tk_1 != null)
				{
					GeoTriangle t = mapping.triangles[tk_1]; // TODO: perhaps find a way to abstract triangleArea to take a triangle
					GeoVertex v0 = mapping.getVertex(t.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
					GeoVertex v1 = mapping.getVertex(t.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
					GeoVertex v2 = mapping.getVertex(t.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
					sharedArea += Geometry.triangleArea(v0, v1, v2);
				}
				if (tk_2 != null)
				{
					GeoTriangle t = mapping.triangles[tk_1];
					GeoVertex v0 = mapping.getVertex(t.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
					GeoVertex v1 = mapping.getVertex(t.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
					GeoVertex v2 = mapping.getVertex(t.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
					sharedArea += Geometry.triangleArea(v0, v1, v2);
				}
				shared.put(nk, sharedArea);
				SWi += sharedArea;
			}
			sharedTriangleAreas.put(vk, shared);
			totalTriangleAreas.put(vk, SWi);
			gv.setNormal(n);
		}
		
		// now we compute the curvature values
		for (int vk = 0; vk < mapping.vertices.length; vk++)
		{
			GeoVertex gv = mapping.getVertex(vk, CURRENT_SMOOTHNESS_LEVEL); // possible smoothed vertex (v_i)
			
			float[] n = gv.geNormal();
			double SWi = totalTriangleAreas.get(vk);
			double[][] tensor = new double[3][3]; // curvature tensor
			for (Integer nk: gv.neighbourkeys)
			{
				GeoVertex gvn = mapping.getVertex(nk, CURRENT_SMOOTHNESS_LEVEL); // v_j
				double[] vij_diff = Geometry.vectorMinus(gv, gvn);
				double[] vji_diff = Geometry.vectorMinus(gvn, gv);
				
				// compute the vector Tij = [(I - N*N^t)(vi - vj)]/|(I - N*N^t)(vi - vj)|
				double[] Tij = new double[3];
				Tij[0] = (1-n[0]*n[0])*vij_diff[0] - n[1]*n[0]*vij_diff[1] - n[2]*n[0]*vij_diff[2];
				Tij[1] = -1*n[0]*n[1]*vij_diff[0] + (1-n[1]*n[1])*vij_diff[1] - n[2]*n[1]*vij_diff[2];
				Tij[2] = -1*n[0]*n[2]*vij_diff[0] - n[1]*n[2]*vij_diff[1] + (1-n[2]*n[2])*vij_diff[2];
				
				double Tij_Mag = Geometry.vectorMagnitude(Tij);
				Tij[0] /= Tij_Mag;
				Tij[1] /= Tij_Mag;
				Tij[2] /= Tij_Mag;
				
				// now compute Kij = [ 2N^t (vj - vi) ] / |(vj - vi)|^2
				double Kij = 2*n[0]*vji_diff[0] + 2*n[1]*vji_diff[1] + 2*n[2]*vji_diff[2];
				double vji_diff_mag = Geometry.vectorMagnitude(vji_diff);
				Kij /= (vji_diff_mag*vji_diff_mag);
				
				// compute M += Wij Kij Tij Tij^t
				double Wij = sharedTriangleAreas.get(vk).get(nk) / SWi;
				double weight = Wij * Kij;
				tensor[0][0] += (Tij[0]*Tij[0] * weight); // row 0
				tensor[0][1] += (Tij[1]*Tij[0] * weight);
				tensor[0][2] += (Tij[2]*Tij[0] * weight);
				
				tensor[1][0] += (Tij[0]*Tij[1] * weight); // row 1
				tensor[1][1] += (Tij[1]*Tij[1] * weight);
				tensor[1][2] += (Tij[2]*Tij[1] * weight);
				
				tensor[2][0] += (Tij[0]*Tij[2] * weight); // row 2
				tensor[2][1] += (Tij[1]*Tij[2] * weight);
				tensor[2][2] += (Tij[2]*Tij[2] * weight);
			}
			
			// now that we have the tensor, we need to compute the householder transformation 
//			double[] W1 = {1-n[0], -n[1], -n[2]};
//			double[] W2 = {1+n[0], n[1], n[2]};
//			
//			double W1_mag = Geometry.vectorMagnitude(W1);
//			double W2_mag = Geometry.vectorMagnitude(W2);
//			
//			int sign = 1;
//			if (W1_mag > W2_mag) sign = -1;
//			
//			double[] W = new double[3];
//			if (sign < 0)
//			{
//				W[0] = W1[0]/W1_mag;
//				W[1] = W1[1]/W1_mag;
//				W[2] = W1[2]/W1_mag;
//			}
//			else
//			{
//				W[0] = W2[0]/W2_mag;
//				W[1] = W2[1]/W2_mag;
//				W[2] = W2[2]/W2_mag;				
//			}
//			
//			double[][] Q = new double[3][3];
//			Q[0][0] = 1 - 2*W[0]*W[0];
//			Q[0][1] = - 2*W[1]*W[0];
//			Q[0][2] = - 2*W[2]*W[0];
//			
//			Q[1][0] = - 2*W[0]*W[1];
//			Q[1][1] = 1 - 2*W[1]*W[1];
//			Q[1][2] = 1 - 2*W[2]*W[1];
//			
//			Q[2][0] = - 2*W[0]*W[2];
//			Q[2][1] = - 2*W[1]*W[2];
//			Q[2][2] = 1 - 2*W[2]*W[2];
//			double[][] Q_t = Matrix_3X3.transpose(Q);
//			double[][] curves = Matrix_3X3.multiply(Q, tensor);
//			curves = Matrix_3X3.multiply(curves, Q_t);
//			double k1 = (3*curves[1][1] - curves[2][2]);
//			double k2 = (3*curves[2][2] + curves[1][1]);
			
			// once we have computed the tensor, extract the principle curvature values (tensor[1][1] and tensor[2][2] == principle curvatures??)
			double k1 = (3*tensor[1][1] - tensor[2][2]);
			double k2 = (3*tensor[2][2] + tensor[1][1]);

			if (Math.abs(k1) < 0.00001) k1 = 0;
			if (Math.abs(k2) < 0.00001) k2 = 0;
			gv.principle_curvatures[0] = k1;
			gv.principle_curvatures[1] = k2;
			gv.curvatureMean = (k1 + k2)/2;
			gv.curvatureGaussian = k1*k2;
			
			pc_gaus_min = Math.min(gv.curvatureGaussian, pc_gaus_min);
			pc_gaus_max = Math.max(gv.curvatureGaussian, pc_gaus_max);
			
			pc_mean_min = Math.min(gv.curvatureGaussian, pc_mean_min);
			pc_mean_max = Math.max(gv.curvatureGaussian, pc_mean_max);
		}
		
		System.out.println(" done! [Gaus min="+pc_gaus_min+" , max="+pc_gaus_max+"]\t[Mean min="+pc_mean_min+" , max="+pc_mean_max+"]"); 
	}
	
	// compute k1 and k2 via steiners method
	public void computePrincipleCurvatures_Steiner()
	{
		System.out.print("Computing Principle Curvatures (Steiner)... ");
		for (GeoVertex v: mapping.vertices)
		{
			GeoVertex gv = mapping.getVertex(v.id, CURRENT_SMOOTHNESS_LEVEL);
			// compute a 1-ring neighbourhood to determine the curvature tensor from
			ArrayList<GeoEdge> edge_set = new ArrayList<GeoEdge>();
			for(Integer nk: v.neighbourkeys)
			{
				GeoEdge e = new GeoEdge(gv, mapping.getVertex(nk, CURRENT_SMOOTHNESS_LEVEL));
				edge_set.add(e);
			}
			
			double[] curvature = Curvatures.computeAnisotropicCurvature(edge_set, mapping, CURRENT_SMOOTHNESS_LEVEL);
			// find kmax as it is the curvature value which represents roughness
			double kmax = Math.max(curvature[0], curvature[1]);
			kmax = Math.max(kmax, curvature[2]);
			
			double kmin = Math.min(curvature[0], curvature[1]);
			kmin = Math.min(kmin, curvature[2]);
			
			gv.principle_curvatures[0] = kmax;
			gv.principle_curvatures[1] = kmin;
			gv.curvatureMean = (kmax + kmin)/2;
			gv.curvatureGaussian = (kmax * kmin);
			
			if (Double.isNaN(gv.curvatureMean)) gv.curvatureMean = 0.0;
			if (Double.isNaN(gv.curvatureGaussian)) gv.curvatureGaussian = 0.0;
		}
		System.out.println(" done!");
	}
	
	
	public void computeNRingCurvatureMap(int n)
	{
		computePrincipleCurvatures_Taubin(); // first we need to have curvature values for all the vertices
		
		System.out.print("Computing "+n+"-Ring curvature average... ");
		
		double maxRingCurvature = Double.MIN_VALUE;
		double minRingCurvature = Double.MAX_VALUE;
		for (GeoVertex v: mapping.vertices)
		{
			GeoVertex gv = mapping.getVertex(v.id, CURRENT_SMOOTHNESS_LEVEL);
			gv.curvatureAvgRing = 0.0;
			
			int number = 0;
			HashSet<Integer> visited = new HashSet<Integer>();
			ArrayList<ArrayList<GeoVertex>> ring_set = findRingNeighbourhood(gv, n, visited);
			for (ArrayList<GeoVertex> ring: ring_set)
			{
				for (GeoVertex r: ring)
				{
					gv.curvatureAvgRing += r.curvatureMean;
					number++;
				}
			}
			gv.curvatureAvgRing /= number;
			
			maxRingCurvature = Math.max(gv.curvatureAvgRing, maxRingCurvature);
			minRingCurvature = Math.min(gv.curvatureAvgRing, minRingCurvature);
		}
		System.out.println("done! [min:"+minRingCurvature+", max="+maxRingCurvature+"]");
	}
	
	
	public void computeNRingSimilarityMap(int n)
	{
		computeNRingCurvatureMap(n); // first compute the n-ring curvature map
	}

	
	// returns the current pointset that is used to render the mesh
	public ArrayList<Point_dt> getCurrentPointSet()
	{
		ArrayList<Point_dt> points = new ArrayList<Point_dt>();
		for(GeoVertex g: mapping.vertices) 
		{
			GeoVertex gv = mapping.getVertex(g.id, CURRENT_SMOOTHNESS_LEVEL);
			double x = gv.x()*MESH_ZOOM;
			double y = gv.y()*MESH_ZOOM;
			double z = gv.z()*MESH_ZOOM;
			
			points.add(new Point_dt(x, y, z));
		}
		
		return points;
	}
	
	/**
	 * Return the elements of the mesh as a list of STLTriangles. Ignore the color attribute.
	 * @return
	 */
	public ArrayList<STLTriangle> getTriangleSet()
	{
		ArrayList<STLTriangle> triangles = new ArrayList<STLTriangle>();
		for (int i = 0; i < mapping.triangles.length; i++)
		{
			GeoTriangle t = mapping.triangles[i];
			GeoVertex v1 = mapping.getVertex(t.vertexkeys[0], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex v2 = mapping.getVertex(t.vertexkeys[1], CURRENT_SMOOTHNESS_LEVEL);
			GeoVertex v3 = mapping.getVertex(t.vertexkeys[2], CURRENT_SMOOTHNESS_LEVEL);
			
			STLTriangle stl_t = new STLTriangle();
			float[] stl_v1 = new float[3];
			float[] stl_v2 = new float[3];
			float[] stl_v3 = new float[3];
			float[] stl_normal = Geometry.computeSurfaceNormal(v1, v2, v3);
			
			stl_v1[0] = (float) v1.x();
			stl_v1[1] = (float) v1.y();
			stl_v1[2] = (float) v1.z();
			
			stl_v2[0] = (float) v2.x();
			stl_v2[1] = (float) v2.y();
			stl_v2[2] = (float) v2.z();
			
			stl_v3[0] = (float) v3.x();
			stl_v3[1] = (float) v3.y();
			stl_v3[2] = (float) v3.z();
			
			stl_t.v1 = stl_v1;
			stl_t.v2 = stl_v2;
			stl_t.v3 = stl_v3;
			stl_t.normal = stl_normal;
			
			triangles.add(stl_t);
		}
		
		return triangles;
	}

	/* 
	 * Triangulation2DFrame functions: used mainly for debugging
	 */ 
	public GeoVertex selectRandomVertex(Random rand)
	{
		int index = rand.nextInt(mapping.vertices.length);
		GeoVertex gv = mapping.vertices[index];
		return gv;
	}
}
